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Abstract 

Background: The aim of this study was to establish an osteosarcoma (OS) associated protein-protein interaction 
network and explore the pathogenesis of osteosarcoma. 

Methods: The gene expression profile GSE9508 was downloaded from the Gene Expression Omnibus database, 
including five samples of non-malignant bone (the control), seven samples for non-metastatic patients (six of which 
were analyzed in duplicate), and 1 1 samples for metastatic patients (10 of which were analyzed in duplicate). 
Differentially expressed genes (DEGs) between osteosarcoma and control samples were identified by packages in R 
with the threshold of |logFC (fold change)| > 1 and false discovery rate < 0.05. Osprey software was used to 
construct the interaction network of DEGs, and genes at protein-protein interaction (PPI) nodes with high degrees 
were identified. The Database for Annotation, Visualization and Integrated Discovery and WebGestalt software were 
then used to perform functional annotation and pathway enrichment analyses for PPI networks, in which P<0.05 
was considered statistically significant. 

Results: Compared to the control samples, the expressions of 42 and 341 genes were altered in non-metastatic OS 
and metastatic OS samples, respectively. A total of 15 significantly enriched functions were obtained with Gene 
Ontology analysis (P< 0.05). The DEGs were classified and significantly enriched in three pathways, including the 
tricarboxylic acid cycle, lysosome and axon guidance. Genes such as HRAS, IDH3A, ATP6apl, ATP6V0D2, SEMA3F and 
SEMA3A were involved in the enriched pathways. 

Conclusions: The hub genes from metastatic OS samples are not only bio-markers of OS, but also help to improve 
therapies for OS. 

Keywords: Osteosarcoma, Differentially expressed genes, Protein-protein interaction network, Pathway enrichment 
analysis 



Background 

Osteosarcoma (OS), an aggressive malignant neoplasm, 
is the most common primary bone cancer in children and 
adolescents [1]. Generally, OS occurs in the long bones of 
the extremities, such as the femur, the tibia and the hu- 
merus [2]. Previous studies have demonstrated that it is 
possible to cure more than 60% of non-metastatic OS with 
neoadjuvant chemotherapy [3]. However, the 5-year sur- 
vival rate of patients with metastases at diagnosis is less 
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than 30% [4,5]. Although the mortality has been decreased 
by about 1.3% per year as a result of improvement in sur- 
gical techniques and chemotherapy, there is no reliable 
means of predicting outcomes at the time of presentation. 
Thus, identifying potential therapeutic targets is very im- 
portant in the treatment of patients with metastatic OS. 

OS is a heterogeneous disease including metastatic and 
non-metastatic complex etiologies. Extensive efforts have 
been made to explore the underlying pathogenesis and 
molecular mechanisms. As reported [6], Ezrin links the 
plasma membrane to the actin cytoskeleton where it plays 
a pivotal role in the metastatic progression of cancers. 
Ezrin immunoreactivity can be detected in the majority of 
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patients with non-metastatic OS of the extremity [6]. How- 
ever, metastatic cells with an efficiently regulated actin 
cytoskeleton may have a greater fluency and flexibility in 
managing energetic needs during times of stress than non- 
metastatic cells. The dysregulation of Ezrin may lead to an 
urgent but failed attempt by cells to correct this metabolic 
defect [7]. Therefore, Ezrin is higher expressed in meta- 
static cells. There are also many other genes that are chan- 
ged in metastatic OS. Clusterin, collagen 4A3, integrin fi4, 
integrin fi2, vascular endothelial growth factor, connective 
tissue growth factor, integrin-aS, cadherin 2, exc ligand 12, 
exc chemokin receptor 4, cell division cycle 5 like, and 
microRNA215 are found to be upregulated in metastatic 
OS, while cadherin 11, MME, KIAA19S9 and S100A6 are 
downregulated in metastatic OS [8-11]. Furthermore, some 
prognosis factors have been demonstrated as potential 
therapeutic targets, such as reversion-inducing-cysteine- 
rich protein with kazal motifs and urokinase plasminogen 
activator [10,11]. The receptor tyrosine kinase-like orphan 
receptor 2 has also have been identified as being important 
in the development of metastatic OS [8,12]. In spite of the 
great progress in the pathogenesis of OS, the molecular 
mechanisms have not been fully elucidated. 

In this paper, we focus on exploring the molecular 
mechanisms of OS with a computational bioinformatics 
analysis of gene expression. We anticipate that candidate 
genes identified by our method may provide a new thera- 
peutic approach for OS. 

Methods 

Affymetrix microarray 

The gene profile of GSE9508 was downloaded from the 
Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/ 
geo/), which was based on the platform of Agilent- Whole 
Human Genome Oligo Microarray G4112A, condensed. 
This expression dataset was deposited by Endo-Munoz 
[13]. Biopsies were collected at the time of diagnosis, 
before preoperative chemotherapy, near the tumor-bone 
interface and were representative of the tumor. A total of 
39 samples were derived from 22 patients (11 males and 
11 females), including five samples of non-malignant bone 
(the control), seven samples for non-metastatic patients 
(six of which were analyzed in duplicate), and 11 samples 
for metastatic patients (10 of which were analyzed in du- 
plicate). Non-malignant bone was collected from patients 
presenting for hip or knee replacement surgery. 

Data processing and identification of differential 
expression analysis 

Firstly, data at the probe-level in CEL format were con- 
verted into expression measures and median normalization 
was performed by R language [14,15]. Then, differentially 
expressed genes (DEGs) between the samples were identi- 
fied with the Limma package [16] in R. The Benjamini & 



Hochberg [17] method was used to adjust the raw P-values 
into the false discovery rate (FDR). |logFC (fold change) | > 1 
and FDR < 0.05 were chosen as the cut-off criteria. 

Construction of protein-protein interaction network 

The interactions between the DEGs in the groups were 
analyzed and the protein-protein interaction (PPI) net- 
works were constructed with Osprey using the data from 
the General Repository for Interaction Datasets [18,19] 
and Bio-molecular Interaction Network Database [20]. A 
total of 159 interaction pairs were analyzed. 

Gene ontology enrichment analysis 

Gene ontology (GO) analysis has been used frequently 
in functional studies of large-scale genomic or transcrip- 
tomic data [21,22]. GO provides three structured networks 
of defined terms, namely biological process, molecular 
function, and cellular compartment [23]. To functionally 
classify PPI nodes with a high degree (number of inter- 
actions), we performed biological process ontology. We 
performed GO enrichment analysis with The Database 
for Annotation, Visualization and Integrated Discovery 
(DAVID), which provides functional interpretation for 
researchers to study biological information behind large 
lists of genes [24]. A P- value less than 0.05 was consid- 
ered statistically significant. 

Pathway enrichment analysis 

We utilized the WebGestalt web site (web-based gene set 
analysis toolkit, http://genereg.ornl.gov/webgestalt) [25,26] 
for pathway enrichment analysis. A P-value less than 0.05 
was considered as the cut-off criterion. 

Results 

Differential gene expression in osteosarcoma samples 

After data normalization (Figure 1), a total of 42 and 
341 DEGs with the FDR < 0.05 and |logFC| > 1 were iden- 
tified between non-metastatic or metastatic samples and 
the control. These findings indicated that genetic abnor- 
malities in metastatic samples were significantly higher 
than in non-metastatic samples. 

Protein-protein interaction network 

Interaction networks of DEGs were constructed with 
Osprey (Figure 2A). As a result, interaction relationships 
of 145 DEGs were obtained. Several PPI nodes were 
high in degrees in the network, such as POLE (10), 
HRAS (9), H3F3, HSPA4, ADCY2 (8), IDH3A, ATP6apl 
and ATP6V0D2. All these nodes were obtained from 
metastatic samples. A high level of connectivity is shown 
in the network (Figure 2A). All of the five key hub por- 
teins (POLE (10), HRAS (9), H3F3, HSPA4, ADCY2 (8)) 
and connections were then deleted. As shown in Figure 2B, 
integrity of the network was disrupted due to the absence 
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Figure 1 Normalization of gene expressions in metastatic, non-metastatic and control samples. The medians (black lines) are almost at 
the same level, indicating a good level of standardization. 



of key hubs and connections between the nodes. Com- 
pared to the previous network, isolated node clusters were 
increased in the new network. 

Gene ontology analysis 

We performed functional enrichment analysis using the 
online biological classification tool DAVID, and obtained 
15 significantly enriched functions (Table 1, P< 0.05). The 
PPI nodes of high degrees were classified with GO analysis 
according to their functions. The most remarkable func- 
tion was regulation of adaptive immune response based 
on somatic recombination of immune receptors built from 
immunoglobulin superfamily domains (GO:0002822, P = 
0.01166). The other significant functions included mRNA 
metabolic process (GO:0016071, P = 0.01427), negative 
regulation of molecular function (GO:0044092, P = 
0.02485), nuclear-transcribed mRNA catabolic process 
(GO:0000184, P = 0.02712) and so on. Genes such as C3, 
TBX21, HNRNPL, NCBP2, UPF3A, HRAS, SEMA3F and 
SEM3A were involved in this GO analysis. 

Pathway enrichment analysis of the protein-protein inter- 
action network 

A total of three significant pathways with P-value less 
than 0.05 were enriched (Table 2). The most significant 



pathway was the citrate cycle (hsa00020, tricarboxylic acid 
cycle (TCA) cycle) with a P-value of 0.006557. Three 
genes were involved in the TCA cycle, including SUCLA2, 
OGDH and IDH3A. The other two significant pathways 
were the lysosome (P = 0.007354) and axon guidance {P = 
0.009702). Genes associated with the lysosome pathway 
were IDS, ATP6AP1, GALNS, ATP6V0D2 and AP3B1. Five 
genes were involved in axon guidance, including HRAS, 
RGS3, SEMA3F, SEMA3A, APHB1 and EPHB1. 

Discussion 

In the present study, we investigated gene expression 
profile in OS patients and explored the potential mo- 
lecular mechanisms using computational bioinformatics 
methods. Results showed that expressions of 42 and 341 
genes were altered in non-metastatic OS and metastatic 
OS samples, respectively. The DEGs were significantly 
enriched in three pathways, including the TCA cycle, 
lysosome and axon guidance. PPI nodes such as HRAS, 
IDH3A, ATP6apl and ATP6V0D2 were of high degrees 
in the network, all of which were obtained from meta- 
static samples. 

IDH3A, a mitochondrial NAD-dependent isocitrate 
dehydrogenase involved in the TCA cycle, was identified 
to exert an effect on OS in our study. OS cells have been 
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Figure 2 Interaction networks of DEGs constructed with Osprey. (A) Interaction networks of differentially expressed genes constructed with 
Osprey. (B) The integrity of the network was disrupted due to the absence of five key hub proteins (POLE (10), HRAS (9), H3F3, HSPA4, ADCY2 (8)) 
and connections between the nodes. 



demonstrated to have an impaired TCA cycle and par- 
tially compensate for the mitochondrial defect [27]. This 
gene was differentially expressed in metastatic samples 
compared with the non-metastatic and non-malignant 
samples, suggesting that IDH3A is closely related to mi- 
gration of metastatic OS cells. In the present study, ex- 
pression of oncogenic HRAS was increased in metastatic 



samples. This result was in accordance with a previous 
study where expression of oncogenic HRAS induced per- 
manent growth arrest in human OS cell lines [28]. 

Several genes which are specifically found in metastatic 
OS cells are associated with inhibition of osteoclast differ- 
entiation in osteosartic patients, including ATP6apl (Ac45), 
ATP6V0D2 (an isoform of V-ATPase subunit d) and SEMA 
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Table 1 Functional enrichment of differentially expressed genes between metastatic samples, non-metastatic samples 
and non-malignant samples (P<0.05) 



Category 


Term 


Count 


P-value 


Genes 


GO:0002822 


Regulation of adaptive immune response based on somatic 
recombination of immune receptors built from immunoglobulin 
superfamily domains 


4 


0.01166 


C3, TBX21, RARA, FOXP3 


GO:0002819 


Regulation of adaptive immune response 


4 


0.01225 


C3, TBX2 1, RARA, FOXP3 


GO:0016071 


mRNA metabolic process 


9 


0.01427 


HNRNPL, NCBP2, UPF3A, C170RF71, RNMT, 
GTF2F1, 'sYNCRIP, SNRPC, SF3A 1 


GO:0044092 


Negative regulation of molecular function 


8 


0.02485 


MDFI, KDM 1A, ADCY2, RGS3, WWP2, BUB1B, 
FOXP3, ANGPTL4 


GO:0000184 


Nuclear-transcribed mRNA catabolic process, nonsense- 
mediated decay 


3 


0.02712 


NCBP2, UPF3A, C170RF71 


GO:0007626 


Locomotory behavior 


7 


0.03055 


HRAS, CCR3, SEMA3F, MY015A, CHRNA4, 
SFMA3A, ESPN 


GO:0051129 


Negative regulation of cellular component organization 


5 


0.03381 


SEMA3F, RDX, CLASP2, SEMA3A, FOXP3 


GO:0048841 


Regulation of axon extension involved in axon guidance 


2 


0.03386 


SEMA3F, SEMA3A 


GO:0048843 


Negative regulation of axon extension involved in axon 
guidance 


2 


0.03386 


SEMA3F, SEMA3A 


GO:0006370 


mRNA capping 


2 


0.03386 


NCBP2, RNMT 


GO:0000956 


Nuclear-transcribed mRNA catabolic process 


3 


0.03609 


NCBP2, UPF3A, C170RF71 


GO:0045621 


Positive regulation of lymphocyte differentiation 


3 


0.03609 


RARA, INPP5D, AP3B1 


GO:0051249 


Regulation of lymphocyte activation 


5 


0.03845 


TBX21, RARA, INPP5D, FOXP3, AP3B1 


GO:0009452 


RNA capping 


2 


0.04215 


NCBP2, RNMT 


GO:0050919 


Negative chemotaxis 


2 


0.04215 


SEMA3F, SEMA3A 



(semaphoring). ATP6apl, an accessory subunit of vacuolar- 
type H + -ATPases (V-ATPases), is identified to be involved 
in the progress of OS formation. Significantly reduced 
normal osteoclast formation and inhibition of osteoclast 
differentiation are found in Ac45 decreased osteoclast. Ac45 
knockdown leads to the decrease in the number of osteo- 
clasts, which does not result from abnormal apoptosis but 
from decreased osteoclast precursor cell proliferation and 
fusion. This decrease is partially due to the downregulation 
of extracellular signal-regulated kinase phosphorylation and 
FBJ OS oncogene (c-fos), nuclear factor of activated T-cells, 
cytoplasmic 1 (NFATcl), and Transmembrane 7 superfam- 
ily member 4 expression. Currently, ATP6apl is regarded 
as a novel therapeutic target for osteolytic disease [29]. 



Table 2 Pathway enrichment in the protein-protein inter- 
action network using the WebGestalt 



Category 


Count 


P-value 


Differentially expressed 
genes 


hsa00020: Citrate cycle 
aCA cycle) 


3 


0.006557 


SUCLA2, OGDH, IDH3A 


hsa04142: Lysosome 


5 


0.007354 


IDS, ATP6AP1, GALNS, 
ATP6V0D2, AP3B1 


hsa04360: Axon 
guidance 


5 


0.009702 


HRAS, RGS3, SEMA3F, 
SEMA3A, EPHB1 



A total of three significant pathways with P< 0.05 were enriched. TCA, 
tricarboxylic acid cycle. 



ATP6V0D2 was demonstrated to downregulate the mi- 
gration of metastatic OS cells in this study. Indeed, it 
has been identified to mediate osteoclast differentiation 
[30] and bone formation [31,32]. ATP6v0d2, involved in 
resorption of the bone extracellular matrix, is induced 
by receptor activator of NF-kB ligand (RANKL) which is 
necessary for osteoclast formation. Thus, ATP6v0d2, as 
well as NFAT-cl, TRAP, DC-STAMP and cathepsin K, is 
considered one of the osteoclast marker genes [33,34]. 

We also observed the differential expression of sema- 
phoring family genes (SEMA3F, SEMA3A), a family of 
proteins that exert a profound impact on cancers [35]. 
SEMA3F binds to neuropilin-2 (NRP2) and acts as a 
tumor suppressor [36]. NRP2 has been reported to affect 
the increased vascularization and correlate with advanced 
tumor progression and prognosis in OS [35-38]. Poorer 
prognosis is related to OS with NRP2 rather than those 
without NRP2 [38]. In contrast, SEMA-3A binds to NRP1 
which induces the collapse of neuronal growth cones. Fur- 
thermore, SEMA-3A is induced by RANKL, affecting the 
osteoblastic profile of OS. This demonstrates that target- 
ing RANK-positive OS may be a novel therapeutic ap- 
proach [39]. Overall, metastatic cells have been reported 
more migratory on the premise of bone marrow factors 
than non-metastatic cells. Furthermore, it is the inherent 
migratory ability and the sensitivity to osteoclast-mediated 
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inhibition of migration that shows the major differences 
between metastatic and non-metastatic OS cells [13]. 

Although these PPI nodes were isolated and suspected 
to be involved in the mechanism of metastatic OS, many 
of them are not shown to participate in this process in 
experimental research. Further research will be needed 
to investigate the role played by the genes present in the 
mechanism of metastatic OS. 

Conclusions 

We have demonstrated 42 and 341 DEGs in non- 
metastatic and metastatic OS samples compared to the 
control samples. Significantly enriched pathways (in- 
cluding the TCA cycle, lysosome and axon guidance) 
are dysfunctional in metastatic OS cells. However, the 
specific markers at the molecular level in the pathogen- 
esis of OS are not fully understood. More studies are 
warranted to further explore the molecular mechanisms 
of OS and to apply the molecular genetic diagnosis to 
the clinic. 
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